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We present a simple, general purpose, quantum Monte-Carlo algorithm for out-of-equilibrium 
interacting nanoelectronics systems. It allows one to systematically compute the expansion of any 
physical observable (such as current or density) in powers of the electron-electron interaction cou¬ 
pling constant U. It is based on the out-of-equilibrium Keldysh Green’s function formalism in 
real-time and corresponds to evaluating all the Feynman diagrams to a given order U n (up to 
n = 15 in the present work). A key idea is to explicitly sum over the Keldysh indices in order to 
enforce the unitarity of the time evolution. The coefficients of the expansion can easily be obtained 
for long time, stationary regimes, even at zero temperature. We then illustrate our approach with an 
application to the Anderson model, an archetype interacting mesoscopic system. We recover various 
results of the literature such as the spin susceptibility or the ’’Kondo ridge” in the current-voltage 
characteristics. In this case, we found the Monte-Carlo free of the sign problem even at zero temper¬ 
ature, in the stationary regime and in absence of particle-hole symmetry. The main limitation of the 
method is the lack of convergence of the expansion in U for large U, i.e. a mathematical property of 
the model rather than a limitation of the Monte-Carlo algorithm. Standard extrapolation methods 
of divergent series can be used to evaluate the series in the strong correlation regime. 


The field of electronic correlations is largely dominated 
by applications to strongly correlated materials such as 
high-T c superconductors or heavy fermions. As a re¬ 
sult, the large effort made by the community to build 
new numerical techniques to address correlations aims 
chiefly at reaching strongly correlated regimes for sys¬ 
tems whose one-body dynamics is rather simple (the 
archetype of these systems being the Hubbard model). 
There are, however, many situations where the correla¬ 
tions are either small or moderate, yet their interplay 
with one-body dynamics might be very interesting. Ex¬ 
amples include, for instance, the zero-bias anomaly in dis¬ 
ordered svstem^S, the Fermi- edge singularity in a quan¬ 
tum dolP, a Kondo impurity embedded in an electronic 
interferometer 3 and possibly the 0.7 anomaly in a quan¬ 
tum point contact. While for a few situations, e.g. zero¬ 
dimensional (Kondo effects) and one-dimensional (Lut- 
tinger liquids) syste ms there exist exact analytical and 
numerical techniqueJ^, the vast majority of these prob¬ 
lems remains elusive to theoretical approaches. The aim 
of this article is to design a technique that could ad¬ 
dress moderate interactions for a large variety of out-of- 
equilibrium situations. 

A natural route for dealing with electron-electron inter¬ 
actions is to compute the expansion of physical quantities 
in powers of the interaction coupling constant, denoted 
hereafter by U. This expansion is traditionally written 
in terms of Feynman diagrams. One can then compute 
the first orders, or try various resummation strategies 
that have been elaborated in order to choose the relevant 
Feynman diagrams for a given problem.^ From a numer¬ 
ical point of view, systematic expansions in powers of U 
have also been intensively studied. In this context, vari¬ 
ous diagrammatic Monte-Carlo have been developed and 


studiecP^-R which aim at explicitly summing the series 
of Feynman diagrams numerically, for example for the 
self-energy. Concerning quantum impurity models, there 
has been an intense activity in the recent years in the 
development of new continuous (mostly imaginary) time 
quantum Monte-Carlo techniques, based on an expansion 
in U (or around the strong coupling limit). These new 
algorithms are of huge practical value in solving the self- 
consistent impurity problems that arise from the dynam¬ 
ical mean-field theory of correlated bulk system&SHU] 
even though they still suffer from the sign problem. They 
have been extended to the non-equilibrium case in a rel¬ 
atively straightforward way, simply adapting the Monte- 
Carlo method to the Keldysh formalisnflSSSi. However, 
these out-of-equilibrium versions suffer from a severe dy¬ 
namical sign problem, compared to their equilibrium 
counterparts, which has severely limited their usage in 
practice. In particular, they can not reach the long¬ 
time steady-state limit in several regimes of parameters. 
Moreover, the approach of Ref. IT9l and 1201 has only been 
shown to work with sufficient accuracy for an Anderson 
impurity with particle-hole symmetry, i.e. a very special 
point of the phase diagram. More recently, bold diagram¬ 
matic Monte-Carlo for impurity models have also been 
extended to the Keldysh context and used in com bina- 
tion with the master equation for the density matrbP 1 ^ 
to reach longer time. Finally, testing these approaches 
in large systems, even at moderate interaction, remains 
also an open question. 

In this paper, we first present a simple, systematic and 
general Monte-Carlo method to compute the first 10 to 
15 coefficients of the expansion of any physical observ¬ 
able for a nanoelectronic system, in an out-of-equilibrium 
situation. The system can be a nanoscopic system con- 
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nected to leads or a quantum impurity in a (possibly self- 
consistently determined) bath. Our method can be ap¬ 
plied in various non-equilibrium contexts, either at short 
time after a quench, or in the long-time steady-state. In 
particular, it can easily reach the steady-state limit, even 
at zero temperature, as well as any intermediate time. It 
is also not limited to particle-hole symmetric case. The 
software developed can be seen as an extension of the 
Kwant package^ to tackle electron-electron interactions, 
or of the Triqs packag^ 4 to deal with non-equilibrium sit¬ 
uations. Second, we discuss the issue of the summation 
of the perturbative series of the physical quantity, which 
is well-known to be a prominent topic in the quantum 
many-body problem. We will show that in the out-of- 
equilibrium Anderson model, for the parameters studied 
here, the current through the dot or the density on it 
have a finite apparent radius of convergence as a func¬ 
tion of U. We will also show that by simply modify¬ 
ing the quadratic part of the action (i.e. playing with 
the so-called a term irP3), one can significantly extend 
the radius of convergence, hence in practice compute for 
higher values of the interaction. Finally, we show that 
extrapolation technique for divergent series, e.g. Lin- 
delof method, can also significantly improve the range of 
applicability of our method. 


In section [TJ we summarize our method and explain the 
main differences between our work and previous ones. 
This short section is mostly for QMC experts and can 
be skipped for people new to the field. Section [IT] in¬ 
troduces our models and notations as well as the ba¬ 
sic many-body perturbation expression that forms our 
starting point. This expression relates the interacting 
observables (such as current or magnetization) to the 
non-interacting Green’s function of the system. Section 
|III| discusses how to obtain the latter one, a prerequi¬ 
site to any QMC scheme. While this step is relatively 
straightforward for simple impurity problems (the vast 
majority of the systems considered so far), its general¬ 
ization to non-trivial geometries requires some care, or 
can become a computationally prohibitive task. Section 
ED discusses a direct calculation of the first few orders 
of the interaction expansion by a brute force numerical 
integration. The discussion of the structure of the func¬ 
tions to be integrated will lead to a key insight for the 
Monte-Carlo. Section [V] describes our QMC algorithm. 
In Section [VT] we use the QMC algorithm to calculate the 
first 10-15 terms in the expansion in powers of the inter¬ 
action strength of the local charge on an Anderson impu¬ 
rity in equilibrium. We analyze the radius of convergence 
of the series in presence/absence of a mean-field term in 
the non-interacting Hamiltonian. In section [V 11 [ we use 
the method in the ouf-of-equilibrium regime, to obtain 
some results associated with the Kondo effect. The ar¬ 
ticle ends with a discussion and various appendices that 
contain some proofs or technical details. 


I. SUMMARY OF THE APPROACH 

Let us briefly sketch our algorithm and its proper¬ 
ties. Non-QMC experts can skip this part, since its con¬ 
tent will be detailed and explained in the next sections. 
We start with a general Hamiltonian H(t) = Ho(f) + 
I/Hi n t (t) where Ho(f) is a non-interacting quadratic 
Hamiltonian of an infinite system (typically a nano- 
electronic system connected to several electrodes) and 
Hi n t (Q contains the interacting part which is switched 
on at t = 0. We aim at calculating the expansion of an 
observable Q (say the current or the local occupation of 
an orbital) in powers of U: 

+oo 

Q(U) = £ QnU n (1) 

n—0 

The Q n are given by many-body perturbation theory in 
the Keldysh formalism in the form of a multi-dimensional 
integral of a determinant, according to Wick’s theorem: 

Qn = XV( C ")detM„(C n ) (2) 

C n 

The sum contains an n-dimensional integral over in¬ 
ternal times Ui £ [0, t] as well as a sum over n Keldysh in¬ 
dices at £ {0,1} and a sum over the different interaction 
matrix elements. W(C n ) contains interaction matrix ele¬ 
ments. detM n (C„) is the determinant of a matrix built 
up with the non-interacting Green’s function of H 0 (t). 
Our algorithm works as follows: 

(1) We compute directly the Q n . In contrast, the 
imaginary-time or real-time 1 ' 18 23 continuous-time algo¬ 
rithms sample the partition function of the problem Z. In 
the real-time Keldysh formalism however, the partition 
function is Z = 1 by construction and as we shall see, its 
sampling is not well suited for obtaining the Q n . Techni¬ 
cally, the integrand detM„(C ra ) is concentrated around 
times Ui which are close to t while in the sampling of Z 
the Ui are spread over the full interval [0, t]. This first 
step ensures that our technique converges well as t —> oo 
and that this limit can be taken order by order in U. 

(2) In the Keldysh formalism, one typically starts from 
a non interacting system, switches on the interaction, 
let the system evolve for a time t , measures the observ¬ 
able and then evolves back to the non-interacting ini¬ 
tial state. The Keldysh indices cq are reminiscent of 
the two evolutions [from 0 to f (a, = 0) and back to 
0 (dj = 1)]. The time evolution is unitary. To keep 
this unitarity order by order, we choose to sum explicitly 
over the Keldysh indices. Hence our algorithm samples 
directly | )C{ ai } W(C n ) det M n (C ra )|. Indeed, performing 
the sum over the Keldysh indices only with the Monte- 
Carlo Markov chain implies that unitarity is only re¬ 
spected on average (i.e. not for a single configuration), 
and in other words relies on the Monte-Carlo to perform 
massive cancellations. Obviously the explicit sum comes 
at an exponential cost: for each Monte-Carlo move one 
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needs to calculate 2 n terms. However, we shall see that 
the gain in signal to noise ratio more than compensates 
for this additional computational cost. In the problems 
that we have computed so far, the dynamical sign/phase 
problem entirely disappears even for large times. From a 
diagrammatic perspective, the summation over Keldysh 
indices also implements the automatic cancellation of dis¬ 
connected Feynman diagrams. 

(3) We compute the Q n coefficients individually rather 
than Q(U). Indeed, once the Q n are obtained, we can 
analyze the convergence of the series for Q(U), which 
is a separate mathematical problem and has nothing to 
do with the Monte-Carlo or any other technique used 
to obtain the Q n . Moreover, in continuous-time algo¬ 
rithms, the interaction very often takes the form of a 
density-density interaction H; nt = (n-j- — a-f)(n^ — cq.) 
and very special values of a a must be used for the com¬ 
putation not to be plagued by the sign problenPSHH. For 
instance in Ref. [2D] the algorithm fluctuates wildly away 
from ota = 1/2. Here we find that the convergence of 
the series for Q{U) depends strongly on the value of a a . 
Note that this is a property of the perturbation series 
and therefore independent of the procedure used to ob¬ 
tain the Q n coefficients. In practice, the dependence of 
the radius of convergence on the a a parameters can be 
used to access larger values of the interaction: we add 
(to H 0 ) and substract (to Hint) an on-site potential to 
our Hamiltonian such that the full Hamiltonian is un¬ 
changed but the series acquires a larger radius of conver¬ 
gence. This step allows us to tackle systems away from 
the particle-hole symmetry point. It is not linked to the 
QMC technique per se but to the choice of the initial 
quadratic Hamiltonian around which one performs the 
interaction expansion. 


II. MODELS AND BASIC FORMALISM 



T 2 , |i 2 


FIG. 1. Sketch of a typical mesoscopic system: a central 
interacting region (red) is connected to several (semi-infinite) 
non-interacting electrodes (blue) with finite temperatures Ti 
and chemical potentials /r;. 


A. Models 

We consider a general time-dependent model for a con¬ 
fined nanoelectronic system connected to metallic elec¬ 
trodes, following the approach of Ref. [27] A sketch of 
a generic system is given in Fig. |T| The Hamiltonian 
consists of a quadratic term and an electron-electron in¬ 
teracting term, 

H(i) = H 0 (t) + I/Hi n t (t) (3) 

where the parameter U controls the magnitude of the 
interaction. The non-interacting Hamiltonian takes the 
following form, 

H 0 (t) = ^H° J .(t)cfc j (4) 

hj 

where cj ( Cj ) are the usual fermionic creation (annihila¬ 
tion) operators of a one-particle state on the site i. The 
site index i is general and can include different kinds 
of degrees of freedom: space, spin, orbitals. A crucial 
aspect is that the number of “sites” is infinite so that 
the non-interacting system has a well-defined density of 
states (as opposed to a sum of delta functions for a fi¬ 
nite system) while interactions only take place in a finite 
region. Typically, the system will consist of a central 
part connected to semi-infinite periodic non-interacting 
leads. The dynamics of such non-interacting systems is 
well known and mature techniques exist to calculate both 
their stationarjf 2 ^ and time-dependent propertied^. The 
interaction Hamiltonian reads 

Hj n t(t) — ^ j Vjjki (^)Cjc 7 -CfcC; (5) 

ijkl 

In contrast to the non-interacting part, it is confined to a 
finite region. We also supposed that the interaction van¬ 
ishes for negative time and is slowly or abruptly switched 
on at t = 0. A typical system described by Eq. ([3| is a 
quantum dot where electrostatic gates confine the elec¬ 
trons in a small, highly interacting, region while the elec¬ 
trodes have high electronic density, hence weak interac¬ 
tions. 

The techniques described below are rather general and 
will be discussed within the framework of Eq. ([3|. The 
practical calculations however will be performed on the 
following Anderson impurity models. Model A corre¬ 
sponds to one interacting site, “0” inside an infinite one- 
dimensional chain, 

+oo 

Ha = ^ ^2 7 iC- jCT c.j + i :(T + h.c. + e d (n T -1- nj.) 

i =—OO <7 

+ E/0(£)n t n ; - h( n f - nj.) (6) 

where 


Her — Cq^Cqct 


( 7 ) 
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e d is the level on-site energy, h the (Zeeman) magnetic 
field and 6{t) is the Heaviside function so that the inter¬ 
action is switched on at t = 0. The hopping parameter 7 \ 
is equal to unity 7 ; = 1 for all sites except 7 _i = 70 = 7 . 
We apply a bias voltage 14 between the two (Left and 
Right) electrodes which are characterized by their chem¬ 
ical potential /x^ = 14 and g R = 0 and temperature T. 
Model B is very close to model A with additional param¬ 
eters <Xf, Oq, 

+ OO 

H s = E ^2lM,^i+i^ + h.c + e d (A f + nj.) 

i = — OO (7 

~h( n t - nj.) + U6{t) (n t - a t ) ( 114 . - a;) ( 8 ) 

One easily realizes that the two models are in fact 
equivalent in the stationary limit for 07 = oq = a: 
Hs(ed, [/, a) = — Ua,U) + Ua 2 . However they 

have very different large U limit at fixed (small) e d ' 
model A corresponds to the degeneracy point between 0 
and 1 electrons on the impurity (where Coulomb blockade 
is lifted) while model B corresponds to the Hondo regime. 
More importantly, the perturbation series in powers of U 
of the same observable will be different between these 
two models, with different convergence radius for fixed 
e d ■ The a parameters have been introduced in Ref. ua 
to improve the sign problem in imaginary-time Quan¬ 
tum Monte-Carlo. An important energy scale for these 
models is the (non-interacting) tunneling rate from the 
impurity to the reservoirs. It is given byT = T^ + T/{ 
with T L/R = 272^/1 - (g, L/R /2) 2 . 


B. Interaction Expansion 

Our starting point for this work is a formal expansion 
of the out-of-equilibrium (Keldysh) Green’s function in 
powers of electron-electron interactions. This is a stan¬ 
dard stepPS which we briefly sketch to introduce our no¬ 
tations. 


The expansion in powers of U can now be performed, 

+°° (_n)n r 

<%•(*,?) = -*E ~^r un E(- 1 ) E ' a * / du ^ dU2 

n =0 H ' { Qi } d 


Using the interaction representation, one defines 
Cj(f) = Uo(0,t)cjUo(t, 0) where Uo(£',£) is the evolu¬ 
tion operator from t to t' associated with Ho- Introduc¬ 
ing the Keldysh index a = 0,1, one defines the contour 
ordering for pairs t = ( t,a ): (f, 0 ) < (£', 1 ) for all t,t\ 
(t, 0) < (t',0) if t < t' and (f, 1 ) < (t',1) if t > t'. 
The contour ordering operator T c acts on products of 
fermionic operators A, B,C ... labeled by various “con¬ 
tour times” £a = (tA, a a), tc ■ ■ ■ and reorder them 
according to the contour ordering: T c (A(tA)B(tB) = AB 
if Ia > ts and T c (A(£a)H(£b) = —BA if Ia < ts- The 
non-interacting contour Green’s function is defined as 

QtjitJ) = -i{T c Ci(t)c](P)) ( 9 ) 

where Cj(f) is just Cj(t), the Keldysh index serving only to 
define the position of the operator after contour ordering. 
The contour Green’s function has a matrix structure in 
a, a 1 which reads 




gfjW) 


( 10 ) 


wher egfj(t,tf), gfj(t,t'), g^(t,t') and g?j(t,t') are respec¬ 
tively the time ordered, lesser, greater and anti-time or¬ 
dered Green’s functions. Efficient techniques to obtain 
these non-interacting objects for large systems will be 
discussed in the next section. Last, one defines the full 
Green’s function G'^(t. ?) with definitions identical to 

the above except that Uo is replaced by U, the evolu¬ 
tion operator associated to the full Hamiltonian H. The 
fundamental expression for G(j-(f, t') reads 


G?.(t,F) = -i(T c e~ ldm ™‘“t(u) £ .(^ e t( ? )) (n) 

where the integral over u is taken along the Keldysh con¬ 
tour, i.e. increasing u for a = 0 and decreasing for a = 1 . 
Hint(xx) is equal to Hi nt (tx) with the operators Cj,cj re¬ 
placed by Ci(u),Cj(u). 


du n {T c H int (tXi ) Hint {U2 ) • ■ ■ Hint {u n )Ci (t)Cj(t')) (12) 


Of particular interest to us are one-particle observables (say current or electronic density) which can be directly 
expressed in terms of the lesser Green’s function at equal times: 

On = (U( 0 , t)cjcjU(£, 0)) = -iGf^t) (13) 

Note at this stage that the following derivation is presented for the one particle correlator, but can be straightforwardly 
generalized to higher correlators. 

To proceed, one evaluates the average of the (large) products of creation/destruction operators using Wick theorem, 
in a form of the determinant of a (2n + 1 ) x (2 n + 1 ) matrix, 


+ OO 


n—0 


{«»} 


duidu2 ... du , 


E 

iijikih 


Viijikili (^l) ' ' ' ^ ^ Vi n j n k n l n (^n) d.6t 

i n j n k n l n 


(14) 
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where M„ is given by 


/ 9 kl is u l> u l) 


1>“2) 

- 9 klj ( u i,?) 

\ 


9hj i( Ul > Ul ) 

gf 1 i 2 (u 1 ,u 2 ) 

- 9f d (u i,t') 


9k 2 iS U 2>“l) 

5fc 2jl (u2,Wl) 

g^ i2 (u 2 ,u 2 ) 

- 9 C k 2 j(u 2 ,t') 


9k n ii( U n,Ul) 


9k n i 2 (u„,U 2 ) 

- 9 C krlj {Un,t') 


szU s i) 

9l njl {Un,Ui) 

9? ni2 (Un,U 2 ) 

- 9f n j(u„,f) 


\ 9iiSt, u l) 


9iijt,U2) 

- 9ij{t,t') 

) 


(15) 


and the zeroth order term is gfj (t. t r ). 


Eq. (14) might look cumbersome at first sight, yet it is 


a compact expression: provided one knows how to cal¬ 
culate non-interacting Green’s functions (which will be 
taken care of in the next section), Eq. (14) expresses 


the full interacting Keldysh Green’s function, hence the 
physical observables, in terms of integrals and sums of 
determinant of known quantities. All that remains is to 
find a suitable numerical way to perform those integrals 
and sums. For a local interaction present on L sites, 
calculating all contributions to order U n in the station¬ 
ary regime corresponds to a numerical complexity of the 
order of L n t n where the measurement time t has to be 
large enough for the effect of the electron-electron inter¬ 
action to be well established. This can be performed us¬ 
ing standard integration routines for the first few orders 
(In section IV we calculate contributions n = 0,1, 2, 3,4 
for model A) but becomes quickly prohibitive for larger 
values of n. For larger orders, a stochastic sampling of 
the integrals is compulsory. 

To simplify the notations, we introduce the notion of 
configuration C n 

Cn — (tl > jl 5 ^1 > ^1) til ? • * • i fn) jn: kn > > ^n) ( 16 ) 

and note 

^ = / du\du 2 ...du n E] ••• E] (17) 

C n ./0<' Ul <•••<u Il <max(M , ) njikih i n j n k„l n 


Introducing, 


V(C n ) = n VU, 


h / 

” 'T> "P 


p= 1 


(18) 


we get the following compact expression: 

+ OO 


G b( i ' ? ) = E inc/n E(- 1 ) E * ai x 

n =0 {a;} 

E^(C„)detM„(C„,{a i }) (19) 


function is unity which translates into 

+ OO 

0 = E i n U n E (-!) Ei “* E det P »(C»> {a <}) 

n=l {m } C„ 

( 20 ) 

where the 2 n x 2 n matrix P„ is identical to M„ with the 
last row and column deleted. Actually, a much stronger 
statement can be made on P„: for any n > 0 and con¬ 
figuration C n , one has, 

E(" l) E ‘ 0i detP„(C n ,{a < })=0 (21) 


The proof is straightforward and standard: one first lo¬ 
cates the largest time in the configuration C ra , say u n . 
When a n goes from 0 to 1, the ordering of u n with re¬ 
spect to the other times is unchanged ( u n is larger than 
all the times on the upper part of the contour and smaller 
than all those on the lower part of the contour), hence 
the contour Green’s functions are unchanged and the ma¬ 
trix P„ is also unchanged. As a result of the (—1)“" sign 
these two contributions cancel each other. 


III. THE NON-INTERACTING GREEN’S 
FUNCTION 

In order to proceed with evaluating the interaction cor¬ 
rections to observables, the first step is an efficient way 
to calculate the various real-time non-interacting Green’s 
functions of the problem. For a small dot problem or 
a DMFT model, this step is easy. For larger systems, 
this question is more delicate, and it has been studied 
extensively®!! and we briefly summarize the main aspects 
here. Note that in Ref.[23 only quantum transport was of 
interest so that contributions coming from bound states 
could have been omitted. Here however, they will have 
to be taken into account properly. 


A. General method 


where the n\ factor has dropped out due to the ordering of Our starting point for calculating non-interacting 
the Ui. Note that in the Keldysh formalism the partition Green’s functions is an expression that relates them to 
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the (Scattering) wave functions in the systeirPH, 

= *E [ Y~ f<x( E )^ocE{t,i)^* aE (t',j) 

a ^ ^ 

+ iY / f(E n )* n (t,i)K(t',j) ( 22 ) 

n 

Here, a labels the various propagating channels of the 
leads, v F a £(f,i) the scattering state at energy E (in the 
electrode) and f a (E) the corresponding Fermi distribu¬ 
tion function, n labels a bound state of energy E n and 
wave functions ’F n . The greater Green’s function g/(t. t') 
is obtained with an identical expression with the Fermi 
functions f(E) replaced by f(E) — l. Efficient techniques 
for calculating the scattering wave functions $„£(f,i) 
have been designed so that these objects can be obtained 
for large systems ( 10 5 siteiJ^j). The bound states contri¬ 
bution was not considered in Ref. (28] and will be discussed 
below. The actual calculations performed in this article 
were restricted to a stationary non-interacting system, 
where the above expression further simplifies into 

a ^ ^ 

+ * E f(E n )t> n (i)** n {j)e- iE ^-^ (23) 

n 

Here again, the stationary wave functions , F a £;(z) are 
standard objects. They are in fact direct outputs of 
the Kwant software 23 which we use for their calculations. 
Once the lesser and greater Green’s functions are known, 
one completes the 2x2 Keldysh matrix with the standard 
relations 

gJj (t, t') =0{t- t')gfj (t, t') + 0(t' - t)gf- (:t , t') (24) 
9% (t, t) =9(1/- t)g> 3 (: 1 , t') + 9(t- t’)g< (t, t') (25) 

To obtain those Green’s function numerically, i.e. for 
many values of t—t r , one needs to perform the integration 
over the energy E many times. In practice, the station¬ 
ary wave functions are calculated once using Kwant and 
cached. The integration itself is performed using stan¬ 
dard numerical routines. For the single site model A or B, 
the above technique in its full generality can be avoided: 
one can simply compute the Green’s function in energy 
analytically and perform a numerical Fourier transform. 
We have checked explicitely that both techniques provide 
identical non-interacting Green’s functions in this special 
case. 


B. Bound states contribution 


The presence of the electrodes in the system is very 
important physically: it provides the system with a re¬ 
laxation mechanism. Mathematically, the integral in 


Eq. (23) mixes nearby energies so that the resulting 


non-interacting Green’s functions decay (and oscillate) 


at large times. However, in presence of a large enough 
confining energy (far from zero ed parameter in model 
A), true bound states can appear in the system. They 
have energies outside of the electrode bands and there¬ 
fore cannot hybridize with the plane waves of the elec¬ 
trodes. They satisfy the stationary Schrodinger equation 
H°’F ra = E n ty n for the infinite system. Upon integrating 
over the electrode degrees of freedom, they satisfy a sim¬ 
pler (yet non-linear) equation for the interacting region 
only: 

H°’F rl + E (£?„)$„ = E n ^ n (26) 


where E (E) is the retarded self-energy due to the elec¬ 
trode. For a practical calculation, we do as follows: first 
we truncate H° and keep the interacting region plus a 
rather large (yet finite) fraction of the electrodes. We di¬ 
agonalize the corresponding finite matrix and locate the 
eigenvalues that are outside the conducting bands of the 
electrodes. These eigenvalues are used as initial guess 
and we compute the bound states by iteratively solving 
Eq. (26) until convergence. Note that there is an easy 
check to make sure that one uses a complete basis of the 
problem: one must have, 

E/f i^(*)i 2 + Ei^(*)i 2 = 1 ( 2? ) 


which is not verified if some bound states are forgotten. 
Note also that in most of this article, we focus on situa¬ 
tions where there are no bound states in the system. This 
can be easily achieved by using leads which have a larger 
bandwidth than that of the central system, so that any 
bound state that could take place there hybridize with 
the continuum of the lead. 


IV. ANALYSIS OF THE FIRST TERMS OF THE 
PERTURBATIVE SERIES 


Knowing how to get the non-interacting Green’s func¬ 
tions, we are now ready to calculate the perturbation 
series. A first, rather naive, technique would consist in 
calculating the integrals in Eq. (141 using a simple dis¬ 
cretization scheme (Simpson in our case). Only the first 
few orders can be obtained that way, at large computa¬ 
tional cost. Nevertheless, it is rather instructive and also 
serves as a check for the QMC algorithms discussed in 
the next section. We focus on model A and compute the 
local charge Q(U) = (n-j- + nj.) at various orders in U n , 


+oo 

Q{U) = Y J QnU n (28) 

n—0 

Fig. [2] shows the resulting Q n {ed) for n = 0...3. With a 
parallel implementation, results for Q 4 can also be ob¬ 
tained (not shown) at important computational cost and 
<5,5 is prohibitive. All these results will be reproduced 
using the quantum Monte-Carlo sampling with a tiny 
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fraction of the computational time. Note that the sta¬ 
tionary results obtained for Q n do not mean that the se¬ 
ries Eq. p8| ) is convergent, but only that its coefficients 
are well defined. 



FIG. 2. Q n as a function of td for n = 0, 1, 2 and 3. 
The calculations are performed using a direct evaluations of 
the integrals in Eq. ( | 14| ) using the Simpson rule for t = 20, 
7 = 0.5 and T = 0. 


It is very instructive to have a look at the quantity 
which is actually integrated to obtain the Q n . Fig. [3] 
shows the integrand of Q 2 for the 4 values of the pair 
of Keldysh indices ( 01 , 02 ). We see that this integrand 
decays slowly as a function of zq — zt 2 and even more 
slowly as zq or zz 2 get away from the time t where the 
charge is measured. The sign of the integrand changes as 
one changes the Keldysh indices. Fig. [3] should be con¬ 
trasted with Fig. 0] which shows the same integrand but 
now summed over the four Keldysh indices. The inte¬ 
grand shown in Fig. [4] now decays fast as zzi or zq gets 
away from t. This observation can be proven and gener¬ 
alized for higher orders: the integrand decays to 0 when 
a group of Ui is far from the time t where the physical 
observable is measured, Cf. Appendix |B| Finally, Fig. [5] 
shows the same as Fig. [3] but for the matrix P 2 associ¬ 
ated with the partition function. Note that for P 2 , the 
sum on the Keldysh indices simply vanishes, so there is 
no analogous Figure as Fig.[4]for P 2 . In the next section, 
we will use these observations to design a better sampling 
strategy for the Monte-Carlo method. 



FIG. 3. Colorplot of the integrand of Q 2 as a function of 
the two times zq and U2 for model A with /zz, = fiR = 0, 
e d = 0, T = 0 and t = 10. The four panels correspond 
to the 4 possible values of the two Keldysh indices an and 
a 2 . The explicit form of the integrand is /(zq 5 U2, zq, a 2 ) = 
-5m(-l)^ ai det M 2 (zq, zi 2 , eq, a 2 ). 



1 0.064 

0.056 

0.048 

0.040 
S 0.032 
- 0.024 

- 0.016 

- 0.008 
— 10.000 


FIG. 4. Same parameters as in Fig. [3] but the integrand has 
now been summed over Keldysh indices. The colorplot repre¬ 
sents /(zq,u 2 ) = *Ea 1 ,a 2 ( _1 ) Ei “ i detM 2 (wi,M 2 ,ai,a 2 ) (/ 
is real). Note that the integrand is now real, positive and 
concentrated around zq = u 2 = t. 


in the Fock configuration space (i.e. that not only sam¬ 
ples the integrals themselves but also samples the various 
orders n within one process). 


V. QUANTUM MONTE-CARLO 

The direct method of the previous section works in 
principle but is limited in practice to very small orders 
due to its prohibitive computational cost. Stochastic 
methods, such as the Metropolis algorithm, can be ex¬ 
tremely efficient at calculating integrals in high dimen¬ 
sions. In this section, we propose a new route to sample 
the interacting series by constructing a Markov process 


A. Sampling strategy 

Our algorithm is inspired by the conclusion of the pre¬ 
vious section. It consists in i) sampling directly the phys¬ 
ical quantity to be computed (and not the partition func¬ 
tion, which is Z = 1 anyway in the Keldysh formalism), 
and ii) summing explicitly over the Keldysh indices to re¬ 
store unitarity (the symmetry between the two Keldysh 

















FIG. 5. Same parameters as in Fig. [3] but the integrand now 
uses the matrix P 2 instead of M 2 , i.e. is associated to the par¬ 
tition function instead of an observable. The colorplot repre¬ 
sents f(ui,U 2 ,ai, < 22 ) = ai det P 2 (iii, M 2 , ai, 02 ) 


contours) for all configurations. Indeed, it is clear from 
the contrast observed between Fig. [4] and Fig. [5] that it 
is a much better choice, since the integration region is in 
the first case well localized around the time t at which the 
quantity is computed. Sampling P 2 would result in sam¬ 
pling large regions of the Fock space which are irrelevant 
to the actual observable. 

We introduce: 


that samples this distribution. A key aspect of the ap¬ 
proach is that O n is sampled for several (at least two) 
values of n simultaneously: in the average <C • • • ^>, n 
varies between n min and n max . Introducing 

c n = Y J W[Cn]\ (33) 

C n 


we find that the probability p n to be in the order n (in 
practice the fraction of the Monte-Carlo spent in config¬ 
urations at order n) is 


Pn — C-n ^qmc/^i 


qmc/ qmc 


(34) 


The normalization of the total probability p n = 1 
provides the partition function in terms of the c n , Z qmc = 
CnUqmc- Note that the c n are by definition indepen¬ 
dent of the QMC technique used to calculate them and 
in particular of the value of U qmc . Co is simply given by 
the non-interacting value of the observable: cq = |g^(0)|. 
The last item to introduce is the probability q n for the 
fluctuating sign in Eq. (32 ) to be +1 (and (1 — q n ) to be 
— 1). Note that M[C p ]/]M[C p ]| = ±1 is always real (Cf 
Appendix [C]) so that we only average fluctuating signs, 
not phases. We note q n = (1 + s n )/2 so that s n is the 
average sign at a given order. s n and p n are the di¬ 
rect outputs of the computations. With these notations, 
•C S pn M[Cp\/ |M[C P ]| ^>= s n p n and one finally arrives at 


On — Cn Sr, 


(35) 


M[C n ] = - i n+1 ^(-l)^ Q W(C n )detM„(C n ,K}) 

{ai} 

(29) 

which is a real number , as proven in Appendix [Cj We 
construct a Markov process that samples the following 
density of probability, 

P[Cn] = ^—\U^ c M[Cn}\ (30) 

^qmc 

where the denominator Z qmc ensures that the probabil¬ 
ity is normalized. We use the notation Z qmc and U qmc to 
show explicitly that the ’’partition function” Z qmc and 
interaction parameter U qmc belong to the QMC tech¬ 
nique. In particular, the physical value of the interaction 
U can (and will) be distinct from the one(s) U qmc used 
in the Monte-Carlo; also the physical partition function 
is Z = 1 in the Keldysh formalism. 

Introducing O n , the contribution to order U n to the 
observable O (see Eq. ©): 

o(y) = Y,°n un ( 31 ) 

n 

the terms of the perturbative expansion are given by 

On/Z qmc =« 6 pn » (32) 

where -C • ■ • stands for the average over the proba¬ 
bility P[C n ]- All is left is to construct a Markov process 


and 


Cn+l 1 Pn+1 

Cn Cq nic p n 


(36) 


which relates the observable ( O n ) to the output of the 
computations ( p n ,s n )■ As always, Monte-Carlo com¬ 
putations only provide ratios between quantities, see 
Eq. (36). Here, we use our knowledge of the non¬ 
interacting observable (hence of Co) to obtain the c n and 
finally the O n recursively. Eq. ( [36| needs to be applied 
for all n up to the maximal order needed, but its evalu¬ 
ation for various n needs not to be done within the same 
QMC run. 


B. Moves and detailed balance 


Before we can actually perform calculations, we are left 
with a last task: designing a random walk that actually 
samples Eq. (30) with n varying between n m ; n and n max . 


This step is very similar to the construction of other 
continuous-time QMC. We introduce two sorts of moves: 
the moves where we increase n by one unit by adding one 
vertex and the moves where one vertex is deleted. The 
algorithm to obtain the configuration C n {i + 1) at step 
* + 1 from the configuration C n (i) goes as follows: 

(i) Move selection. We choose the move n —> n + 1 
(n —> n — 1) with probability p -j- ( pi ) with p\ + P} = 1. 
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When n = ?i m in (n max ) we havep-j- = 1 (p± = 1) otherwise 
we typically choose pf = Pi = 1 / 2 . 

(ii) Move n —> n + 1. We select u n+ 1 uniformly in 
[0 ,t\. We select (i n+ i,j n+ i,k n+1 ,l n+1 ) uniformly among 
the Ny different terms Vjjki- The overall probability to 
propose the move is 

W t du n+1 = ( 37 ) 


i.e. one samples the matrix P„ (instead of M„) and one 
samples the Keldysh indices (instead of summing on them 
exactly and explicitely). In this scheme, an observable O 
(say the charge Q , possibly resolved in spin) is given by 

0(C/ qmc ) = Z qmc « (-1) E * “* » (43) 

while the partition function Z = 1 is given by 


(iii) Move n n — 1. We select the vertex to be 
removed uniformly between [1... n\. The overall proba¬ 
bility to propose the move is 

= ^ (38) 

n 

(iv) Detailed balance. We choose the acceptance 
probability 94 ( 71 ) and 94 ( 71 ) so that it satisfies the detailed 
balance equation 


Wf(n)du„ +1 q^(n)P[C n ]Y[dui = 

i= 1 

W\.(n + 1 ) 94(71 + 1)-P[C n+l] n dui (39) 


n+l 


i= 1 


Using the Metropolis algorithm this leads to 
'Wi(n+ l)P[C n+1 ] 


94 ( 71 ) = min ^ 
94 ( 77 ) = min 


W t (n)P[C n \ 

W^n-^PjC^] 

W 4 («+[£„] 


(40) 

(41) 


Note that, as usu al in continuous-time quantum Monte- 
Carlo method;!^ ! 1 4 U I ( the factor du n+ 1 is present on 
both sides of the detailed balance equation and even¬ 
tually drops, so that the Monte-Carlo can be performed 
directly in the (time) continuum. 


C. Remarks 


l = +mc<(-l ) E ‘ 01 


det P„ 
I det P, 


» 


(44) 


Constructing a Markov process that samples Eq. ( |42 
can be done similarly to the construction presented in 
the previous subsection. The observable 0(t/ qmc ) can be 
estimated by taking the ratio of the above two equations. 

Although this approach has shown some success, one 
can make the following remarks. 

(1). The sign of detP n can fluctuate strongly so 
that the statistical average in Eq. (44) is very small and 
1 /Zqmc suffers from a very bad signal to noise ratio. This 
is the sign problem that plagues Quantum Monte-Carlo 
techniques for fermions. Eq. (21 1 indicates that this prob¬ 


lem is most probably worse in presence of the Keldysh 
indices where one expects wildly fluctuating signs. This 
is shown by the data in Fig. [5j 

(2) . It is not guaranteed that the most probable con¬ 
figurations sampled by Eq. ( fufj ) are also the ones that 
contribute most to 0(U qmc ). On the contrary: the de¬ 
terminant of P„ depends only on the relative positions of 
the Ui with respect to the others (it is essentially a sum of 
terms of the form gf j (u 1 -u 3 )g^ l (u 2 ~u 3 ) ... g^ q (u„-u 6 )) 
and not at all of t. The integrals contributing to 0(U qmc ) 
on the other hand decay when the iq get away from t. 
Hence for large times, the above scheme samples values 
of the Ui very far from t which therefore contribute very 
little to the actual observable. This was shown explicitly 
in Fig. 0 Fig. |4] and Fig. [5] 

(3) . The signal to noise ratio usually deteriorates 
rapidly with t in these algorithms making it difficult to 
reach the stationary regime. 


We now have a complete practical scheme for calculat¬ 
ing many-body perturbation to a given observable. Be¬ 
fore we embark in concrete examples, let us make a few 
remarks. 


1. Comparison with the sampling of the partition function 

Our sampling strategy should be contrasted with the 
usual approach, used for instance in Ref. [501 which has 
its origin in the imaginary-time techniques and where the 
(density of) probability P[C n , {(+] to be in the configu¬ 
ration C n with the Keldysh indices {a^} is given by 

P[C n ,{ai}} = —— \U^ mc V(C n )detP n (C n ,{a 7 })\ (42) 

^qmc 


2. Role of the explicit sum over the Keldysh indices 

The sampling of Z = 1 discussed above does not pre¬ 
serve the symmetry between the two parts of the Keldysh 
contour for a given configuration (it preserves it in aver¬ 
age, as it should): in the Keldysh formalism, one starts 
from a non-interacting density matrix, switches on the in¬ 
teraction for some time t, measures the observable, then 
unwinds the effect of the interaction until one is back to 
the original density matrix. Here however, a given config¬ 
uration might have a few Keldysh indices on one branch 
and the rest on the other, meaning that a typical configu¬ 
ration does not enforce the symmetry of Keldysh indices, 
which reflects unitarity. From a different perspective, 
the corresponding expansion includes all the Feynman 
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diagrams, including the disconnected diagrams although 
they have a vanishing contribution to the observables. 

In our scheme in contrast, we explicitly sum over all 
Keldysh indices. Eq. (21) indicates that all contributions 
from disconnected diagrams explicitly drop from the cal¬ 
culation. We calculate only connected diagrams which 
should be advantageous. From a technical perspective, 
one finds that the quantity M[C n ] is always real (not 
complex, see Appendix [C]) so that one averages signs in¬ 
stead of phases. One expects that the resulting potential 
sign problem is milder than the so-called phase problem 
which originates from averaging a random phase. 

However, our scheme has an obvious drawback: one 
evaluation of M[C n \ corresponds to the calculation of 2" 
determinants, so that this algorithm complexity now in¬ 
creases exponentially with the maximum order n. We 
show below results for up to n = 15 and it is reasonable 
to expect that one could calculate up to n = 20. Usual 
algorithms to calculate determinants have complexities 
which scale as n 3 . We show in appendix [d] that the fast 
updates of the determinants (with complexities that scale 
as n 2 ) can easily be extended to the present case using 
Gray code, so that the overall complexity of our algo¬ 
rithm scales as 2 n n 2 . Actually, ’’mirror” Keldysh con¬ 
figurations have equal contributions (see Appendix [C]) so 
that only 2 n ~ 1 configurations need to be included. Over¬ 
all, we shall see that the additional computational com¬ 
plexity is more than compensated by the important gain 
in signal to noise ratio. 


3. Statistical errors and the sign problem 

In practice, we calculate the moments O n sequentially 
with a separate QMC computation for each value of n 
(typically with n min = n - 1, n max = n or n min = 0, 
Umax = n) so that one gets a fully controlled (optimally 
flat) histogram of the various orders. Starting from Cq 
(known without error), we iteratively compute c n from 
different runs that use different interaction strengths U n , 

^= An = Pn{Un) (45) 

Cn— 1 C'nPn—l\C'n) 

One can use for instance U n such that p n = p n -i = 1/2 
(this is always possible but not strictly necessary as long 
as p n /p n -i remains of order unity). 

Note that a naive scheme where one would try to eval¬ 
uate all the values of c ra in a single run would run into an 
artificial difficulty: in a single run, the histogram is usu¬ 
ally sharply peaked around one value of n and one cannot 
get a good statistics both for n = 1 and n = n max . A 
small statistics in, say n = 1 leads to a large error in the 
estimate of pi which further corrupts the evaluation of 
all c n , hence O n . 

Coming back to our scheme, the error made on the 
estimate of A n is bounded by SA n /A n < 2/y/N# where 
N# is the number of independent points. Hence we find 


that the (one standard deviation) error on c„ is bounded 

by 


Sc,i 2 n 

c n \/N# 


(46) 


which can be controlled to arbitrary precision provided 
n <C y/N#. Hence, the calculation of the c ni which in¬ 
volves only positive numbers, can be done with extremely 
good accuracy. 

The limitation to the precision - the sign problem - 
takes its origin in the average s n of the sign contained 
in Eq. (l32|). Note that in contrast to other techniques, 


this sign is here in the numerator and is not present in 
the denominator, i.e. a small sign does not necessarily 
mean a sign problem: it can simply mean a small value 
of O n . This analysis is close to e.g. RefsIMTTl The error 
made on the sign s n is given by a Bernoulli law (+1 with 
probability q n , — 1 with probability l — q n ), hence is given 
by 2y/[q n (l — q n )]/N^ which is always smaller than 


Ss n < 




(47) 


(the upper bound is reached when the sign s n becomes 
small). Putting everything together, O n = c n s n implies 
that \SO n /O n \ = <5s n /|s n | + Sc n /c n and we arrive at the 
relative error, 


SO n 1 | 2 n 

or -|U7^# + v^# 


(48) 


In this form, it seems that the smaller the sign, the larger 
the error, hence the sign problem. However, one must 
remember that while c n and s n depend on the QMC al¬ 
gorithm, their product c n s n = O n does not, so that the 
error can be recast into 


SOn 

O n 


< 


\On\y/N^ 


2 n 

7^ 


(49) 


In this second form, it becomes apparent that a bad sign 
problem (small s n ) is equivalent to a bad sampling choice 
which leads to a large c n . The behaviour of the error is 
therefore intimately linked with the growth of c n with 
n which itself depends strongly on the actual probabil¬ 
ity sampled. For instance, if one samples the sum over 
Keldysh indices (instead of summing explicitely over the 
indices as we do), one gets identical O n but much larger 
c n and consequently much smaller s n . In fact the c ra 
would increase by more than a factor 2 " so that the 
overall method would be much less efficient than the one 
proposed here. The global relative error therefore con¬ 
tains three contributions, which reflect respectively the 
total computing time (1 /y/N#), the intrinsic physics of 
the problem ( O n ) and the quality of the choice of the 
sampling procedure (c n ). From the above analysis, the 
colorplots shown in Fig. @ Fig. [4] and Fig. [5]take a differ¬ 
ent meaning. Indeed, one can see that C 2 (the integral of 
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the function displayed in Fig. [4]) is rather small: the func¬ 
tion decays rather quickly when iq, U 2 get away from t. If 
on the other hand we have chosen to sample the Keldysh 
indices in addition (as in the standard schemes), then C 2 
would have been the sum of the integral of (the absolute 
value of) the different panels of Fig. [3] One immediatly 
realizes that the signal to noise ratio would have been 
much smaller. 


4- Convergence of the interacting series 

Our approach separates the calculation of the O n from 
the study of (the convergence of) the series 0 (U) = 
J2 n O n U n itself. This could be used to obtain the full 
U dependence of the observable, but more importantly, 
it allows one to disentangle physical aspects (for instance 
the convergence or lack of of the series) from technical 
ones (the calculation of its elements). The convergence of 
the series will be discussed next. In principle, one could 
use various resummation procedures such as Pade ap- 
proximant, Lindelof analytical continuation and/or Borel 
resummation in order to extrapolate the series from its 
first coefficients. An example of such a procedure is given 
in the appendices. 

It is important to note already at this stage that the 
parameter a , as introduced in model B, plays a crucial 
role in the algorithms sampling the partition function at 
equilibrium 12 ! a s only special values of a are free of the 
sign problem. It is also known that the typical pertur¬ 
bation order in those algorithms is strongly reduced by 
using the best value of a. We shall find in the next sec¬ 
tion that the parameter a has a drastic influence on the 
convergence of the series J2 n O n U n but not on the actual 
calculation of the O n itself. 



FIG. 6 . QMC results for model A at 7 = 1/2, T = 0 and 
t = 10 for td, = 0 (green squares) and td = —0.5 (blue circles), 
as a function of the order n. Red diamonds: model B with 
7 = 1/2, T = 0 and t = 10 for td = —0.5 and a = 0.5. Top 
panel: Q n , central panel: c n and bottom panel: average sign 

Sn • 



l/N 


VI. FIRST RESULTS: ANALYSIS OF THE 
SERIES CONVERGENCE 

A. Bare results 


FIG. 7. QMC results for model A at 7 = 1/2, T = 0 and 
t = 10. Top panel: charge Q(N, U ) as a function of U, for 
different N and for td = 0: N = 5 (Black), N = 7 (red), N = 
13 (blue) and N — 15 (green). Bottom panel: Q(N,U ) as a 
function of 1/A for different U: U = 0.2 (Orange), U = 0.4 
(red), U = 0.6 (black) and U = 0.8 (blue). 


As a first application, we compute the interaction cor¬ 
rections to the charge Q on the impurity in model A. 
Fig. [ 6 ] shows the Q n for n up to n = 15 as well as the 
corresponding c n and s n . We find that the magnitudes 
of the Q n do not appear to decrease with n but rather re¬ 
main of order unity \Q n ss 1. This implies (Hadamard’s 
theorem) that the series has an apparent convergence ra¬ 
dius of order unity (apparent because it relies on extrap¬ 
olating the behaviour of the first known coefficients to 
large orders). We can already notice that the curve for 
model B (a = 1/2) has a very different behaviour with 
a fastly decreasing c n hence Q n ; this aspect will be dis¬ 
cussed later in the text. Fig. [7] shows the truncated 


series 

JV-l 

Q(N,U)=J2QnU n (50) 

n =0 

as a function of U (upper panel) and l/N (lower panel). 
We find a nice convergence for U < U* « 0.6 but a diver¬ 
gence beyond, i.e. we cannot access the physics beyond 
U = 0.6 (U ss r/2) by simply summing the series. This 
divergence has nothing to do with the QMC technique 
which is just a way to calculate the Q n - it belongs to 
the physics of the problem. In the following, we will dis¬ 
cuss two ways to bypass this problem: performing the 
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interaction expansion from a different starting point and 
resumming the series by moving its singularities away 
from the expansion point. 


B. Hartree-Fock series 


Before going on with the QMC results, let us analyse a 
simpler, approximate, series which is obtained from the 
Hartree-Fock approximation (which reduces to Hartree 
for model A). This series will be useful in identifying a 
possible source for the observed apparent radius of con¬ 
vergence. Introducing, 


9 R (E) = - 


E — e d — 7 2 {E + iy/4 - E 2 ) 


(51) 


the Fourier transform of the retarded non-interacting 
Green’s function g R 0 (t), the non-interacting charge at 
equilibrium and T = 0 takes the form 

1 I" 0 

Qo = ~ dE Im g R {E) (52) 

7r J — 2 


In the Hartree approximation, one replaces the on¬ 
site energy by its mean-field value. The fully self- 
consistent Hartree would be defined as Q HF (e d ,U) = 
Qo[e d + UQ HF {td,U)\. Here however, we restrict our¬ 
selves to summing the ladder of tadpole diagrams which 
is sufficient to illustrate our point: Q HF (e d ,U ) = 
Qo[e d + UQo(ed)]- The corresponding series is given by 
Q hf (U) = £„ Q RF U n with 

Qn F = — f dE Im [g R (E)] n (53) 

n J—2 

Note that the first two moments are the exact ones: 
Qo F = Qo and Q\ F = Q i ■ On the other hand, from 
the above construction Eq.(|5l|), we find that Q HF (U ) has 
branch cuts in the complex plane, given by 

E-e d -UQ 0 (e d )-'r 2 (E±iV4: - E 2 ) = 0 for -2 < E < 0 

(54) 

so that one expects the series in powers of U to have a 
finite radius of convergence given by the branch closest 
to U = 0, i.e. equal to unity (see the inset of Fig. [ 8 ]). In¬ 
deed, Fig. [ 8 ] shows an example of the partial sums which 
diverges around U = 1. This is very reminiscent to what 
we have found for model A. This calculaiton illustrates in 
a simple and tractable approximation that a finite radius 
of convergence is due to the existence of singularities or 
branch cuts in the complex plane. 


C. Using a different non-interacting problem 

In this section, we show that by playing with the a 
parameter, one can greatly enhance the radius of con¬ 
veyance of the series hence access correlated regimes. 



FIG. 8 . Hartree-Fock series of model A as a function of U: 
Exact curve Q HF (U) (thick line) and partial sums Q HF (N, U) 
for N = 10,20,30, e d = 0 and 7 = 1/2. The series has a 
divergence at U = 1. The inset shows the analytical structure 
of Q hf (U) in the complex plane (Re U, Im U ): branch cut of 
Q hf (U) (thick magenta line) and convergence radius of the 
Hartree-Fock series (dashed line). 


Instead of starting the perturbation from the U = 0 
Hamiltonian, one can incorporate the mean-field treat¬ 
ment into the non-interacting Hamiltonian so that only 
the fluctuations of the interaction need to be taken into 
account in the perturbation. In fact, one can even push 
this idea further and add an arbitrary quadratic Hamil¬ 
tonian to Ho(f) and remove it accordingly from H; nt (t). 
The idea is to start with one-body propagators as close 
as possible to the interacting ones, so that the role of the 
perturbation becomes very weak. Within our current im¬ 
plementation we can easily add an on-site potential Se d 
to the non-interacting Hamiltonian and use a = Se d /U 
in the perturbation where the new parameter U is our 
targetted value of the interaction (see the definition of 
model B in Eq. ([ 8 ])). For U = U, one recovers the orig¬ 
inal model A for t > 0, hence the corresponding results 
in the long-time limit. 

The a parameters are exactly the same as the ones 
used in the interaction expansion continuous-time quan¬ 
tum Monte-Carlo introduced by Rubtso^ 12 } in equilib¬ 
rium. In that algorithm, it was shown 12 ! that the sign 
problem strongly depends on the value of a. It also re¬ 
duces the average order of perturbation of these QMC 
methods Here, we will now show that the apparent ra¬ 
dius of convergence of the interaction expansion strongly 
depends on these a parameters, and we will use this to 
our advantage. 

The technique is illustrated in Fig. [9]first for a value of 
U that could be reached with the initial approach (U = 
0.25) and secondly for a value that could not be reached 
{U = 1). We find that this approach works remarkably 
well: not only we can recover the former results, but we 
can also go to regimes that were not accessible before. 
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As a self-consistency check, we find that the results do 
not depend on Std for U = U. Of course a disadvantage 
of this approach is that one must perform a separate 
computation for each value of U needed. 


Q(U) 



FIG. 9. Role of the non-interacting Hamiltonian around 
which is done the expansion; Q(U) at td = 0 with an extra 
potential Std in the non-interacting Hamiltonian and a corre¬ 
sponding compensating term a = 5td/U in the perturbation. 
The green line corresponds to the original series with Sea = 0 
in its regime of convergence. All curves for 7 = 1/2, T = 0 
and t = 10 . 



FIG. 10. Top panel: absolute value of the charges Q n for 
different values of a,Std, at td = 0, 7 = 1/2, T = 0 and t = 10: 
a = 5td = 0 (Black), a = Std = 0.5 (Blue), a = 8t d = 0.6 
(Green), a = Std = 0.7 (Red). Bottom panel: corresponding 
dependence of the charge Q(N, U = 1) as a function of 1/A. 


The faster convergence of the new series can be seen 
in Fig. |10| where we compute the corresponding series Q n 
versus n (upper panel) as well as the convergence of the 
partial sum Q(N, U = U) versus 1 /A (lower panel). We 
find that only 2 or 3 orders are sufficient to obtain the 
exact result provided one uses a ” starting point” close 
enough to the final solution. A pragmatic way to per¬ 
form the calculation is therefore to optimize the value 
of 5ed so that the corresponding Q n decreases as fast as 
possible. We collect the final curve Q(U) for model A 
at ed = 0 (note that other values are equally accessi¬ 
ble) in Fig. 11 together with the original expansions. We 
find that Q(U) decreases monotonously from Q( 0 ) = 1/2 
to Q(U » 1) ss 1/4. For U > 2 the model is already 
close to its large U limit and fluctuations of charge are 
small. Fig. [TT] is an important result of this paper and 
establishes that the regime of strong interaction can be 
reached from a rather naive perturbative expansion. We 
note that the regime shown in Fig. m is a quite difficult 
one for the technique. Indeed raising either the temper¬ 
ature and/or a bias voltage will result in Green’s func¬ 
tions that decay rapidly with time (exponentially as op¬ 
posed to the algebric decay found at zero temperature) 
and therefore in smaller c n and better signal to noise ra¬ 
tio in the calculations. Fig. m also includes a separate 
calculation performed with an hybridization QMC tech¬ 
nique in imaginary time (dashed line), obtained with the 
algorithm introduced in Ref. Il3l implemented with the 
TRIQS packag^l. We find a perfect match between the 
two techniques which serve as a validation of our tech¬ 
nique and its implementation. 



FIG. 11. Charge Q{U) for model A with td = 0, 7 = 1/2, 
T = 0 and t = 10. Red symbols correspond the use of a 
parameters (different symbols for the same value of U corre¬ 
spond to different 8 t or a). Thin lines correspond to partial 
sum with the bare technique (a = 0): Q(N = 6 , U) (cyan), 
Q(N = 7,U) (green) and Q(N = 8 ,U) (blue). Dashed line: 
Benchmark calculation performed with hybridization QMC in 
imaginary time. Blue squares: extrapolation using the homo¬ 
graphic technique, see sectiorfVI D| 


D. Resummation technique: moving the 
singularities away 

Let us go back to our initial series for model A (with 
a = 0). Our Hartree-Fock analysis suggests that the 
source of divergence lies in the complex plane (magenta 
line in Fig. [ 8 ]). However, the Anderson model does not 
display a transition (hence no singularity) on the real 
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axis for U, which suggests that one should be able to 
analytically continue our result for U > 1. 

There is a large body of literature dedicated to the 
study of diverging or even asymptotic series, see e.g. Refs 
HU M and S3] and various techniques can be used to 
extract useful physical information from them. In Ap¬ 
pendix [Ej we use the Lindelof method to extrapolate the 
series beyond the radius of convergence, for the charge 
Q(U ), and show that we can easily get the correct re¬ 
sult for e.g. U = 1.5. Alternatively, one could use Borel 
resummation followed by a Pade fit (not shown). The 
Borel series, like the Lindelof method, allows one to go 
beyond the initial radius of convergence, but not much 
beyond. It is also rather impractical to do in a controlled 
way. 

Here, we use an alternative route known as the Eu¬ 
ler transform 34 ! The idea is to perform a meromor- 
phic transformation W(U) that sends the singularities 
away from the expansion point while bringing the region 
of interest (real and positive Z7s) closer to zero (with 
W(0) = 0). As an illustration, we choose the homo¬ 
graphic transform 


W(U) = 


bU 

U-a 


(55) 


The method is performed in two steps. First, one obtains 
the expansion for the inverse U(W) of W{U) (defined 
as U[W(U)\ = U): U{W) = U\W + U 2 W 2 + U 3 W 3 ... 
and one constructs the (truncated) series for Q(W) = 
Q[U(W)\. The series Q(W) = Q n W n has a radius of 
convergence Rw which may be much larger than the ini¬ 
tial series if a is close enough to the singularities of Q{U). 
In a second step, one evaluates Q(W) for W = W(U). 
If Rw > W{U ) the result will be convergent so that 
the figure of merit for this transformation is the ratio 
Rw/W(U). In practice, Rw is obtained by performing a 
simple exponential fit for the Q. n oc R^ (see the inset of 
Fig(l2] for an example). It is also very appealing concep¬ 
tually and controlled by the figure of merit Rw/W(U ) 
(in practice Rw/W(U) ~ 1.5 gives very precise results 
with only 7-8 orders). 

Fig |12| gives an example of the resunnnation for U = 6 
(much beyond the radius of convergence of the initial 
series). We find that the figure of merit reaches a rather 
high value Rw/W(U) ~ 2 which allows one to obtain 
the exact result with only a few orders. The increase 
of the radius of convergence is actually quite dramatic. 
In the lower panel we plot the actual result obtained as 
a function of a (the results do not depend on b). We 
find that, in the region where the figure of merit is high 
enough, one observes a nice plateau at the correct value. 
We have also reported the full extrapolated Q{U) curve 
in Fig. 11 (blue squares) which is in very good agreement 
with our reference calculation. This method has two clear 
advantages over other resummation methods. First, the 
technique is controlled with the figure of merit. Second it 
can, in principle, works even in the strong coupling limit. 
For the case studied here, the point U = oo is mapped 


onto W = b = 6 which is well within our convergence 
radius Rw{a = —0.25) « 12 so that we can compute 
Q{U = oo) « 0.18. 



FIG. 12. Homographic transformation of the model A series 
for Q(U ) (same data as FigJrJ with 6 = 6. Upper panel: Ef¬ 
fective radius Rw(a) (ful black line) and homographic trans¬ 
formation W(a, U = 6) (dashed red line) as a function of 
position of the singularity a. Bottom panel: resummation re¬ 
sult Q(U = 6) as a function of a. The dashed line indicates 
the exact result Q fts 0.22. Inset: Q n as a function of n for 
a = —0.25 (squares) and the initial series (circles). 


VII. MORE RESULTS: CURRENT AND 
MAGNETIZATION 

In the previous section, we have computed the charge 
of the impurity, a quantity that misses the important 
physics associated to spin fluctuations, the Kondo effect. 
The Kondo effect in quantum dots has been extensively 
studied and we refer to Ref. [35] and [36] for an account of 
the literature. Here, we merely aim at illustrating our 
method with calculations of current versus voltage and 
magnetization versus field characteristics. 


A. Current 


In Fig. [l3j we present some results obtained for the 
current-voltage characteristics of an Anderson impurity 
connected to two electrodes. 


Fig. 13 shows the resulting 7(14) characteristics for two 


values of e^: = 0 which corresponds to the particle- 

hole symmetric case where the non-interacting problem is 
at resonance; and = —0.5 which breaks this symmetry 
and is non-resonant. We also indicate the perfectly trans¬ 
mitting limit I = 14 with the orange dashed line. The 
ed = 0 case is already very interesting: at small bias, the 
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transmission probability is unity because we are exactly 
at the resonance frequency of the impurity. When one 
increases the interaction strength, we also expect perfect 
transmission but for a totally different reason: the origi¬ 
nal resonance is now shifted to negative energies (—0.6 in 
this instance) but a new one, the Kondo resonance, starts 
to develop at the Fermi level (VJ, < T K ). In practice, we 
find that the I n are extremely small at small voltages so 
that the current is unaffected by the interaction. It is 
only at higher bias than interaction becomes relevant. It 
is interesting to note that this ’’Kondo ridge” which is 
a notoriously difficult regime appears here to be one of 
the most tractable ones. The off-resonant case td = —0.5 
corresponds to a situation where the non-interacting cur¬ 
rent is much smaller than the interacting one so that the 
perturbation series must build the Kondo resonance. We 
find that it does indeed build it as one recovers perfect 
transmission at low bias. Up to the accuracy of the cal¬ 
culations (the error bars are of the order of the symbol 
sizes), the result shown in Fig. 13 is an exact solution 


of the non-equilibrium Anderson model, in its stationary 
regime, with and without particle-hole symmetry. 


for visibility of the higher moments) as a function of time 
t. At short time, the moments grow typically as /„ oc t n 
which simply reflects the fact that I n is an n-dimensional 
integral. After one or two oscillations, they reach their 
stationary value for roughly Tt > 10 but one notices that 
the higher orders converge significantly slower than the 
lower orders. Simulations for large time are not partic¬ 
ularly harder than for short time except for one small 
difficulty: our insertion move has a flat distribution in 
the interval [0, f] but only times close enough to t actually 
contribute to the moments, so that the acceptance proba¬ 
bility for these moves eventually drops when t increases. 
This issue could be circonvaluated by reweighting the 
proposed moves around t. Fig. [15] and Fig. [16] show re¬ 
spectively the corresponding evolution of the signs s n and 
the weights c n as a function of time t. The signs remain 
rather large (the smallest value is S6 = 0.025 asymptoti¬ 
cally) and are not a limitation for the calculations. The 
decrease of the I n with n essentially comes from the cor¬ 
responding decrease of the c n . These Monte-Carlo com¬ 
putations are essentially free of the sign problem. The 
real limitation of the results presented in this section is a 
physical one: the apparent radius of convergence (in U ) 
of the series appears to be finite, of the order of 1.5 — 2 
for the current. 



FIG. 13. Current-voltage characteristics /(14) for model B 
with a = 0.5, 7 = 0.5 and T = 0. The current is in units of 
energy time e/h. Black line and black circles: particle-hole 
symmetric point td = 0 for U — 0 and U — 1.2 respectively. 
Red dot-dashed line and red squares td = —0.5 for U — 0 
and U = 1.2. Dashed line: perfect transmission I = 14. 
Inset: convergence of the results as a function of 1/A where 
A is the total number of moments included for Vj, = 0.3, td = 
—0.5,1/ = 1.2 (red squares) and Vb = 0 .2, td — 0, U = 1.2 
(black circles) 

Let us now analyze a bit further this calculation and 
present the various orders of the expansion of the 
current in powers of U. It is interesting to study how I n 
converge to the stationary value with time. Fig.|T4|shows 
the first seven moments I n (rescaled by U n with U = 3 


FIG. 14. First 6 moments of the currents I„U n versus time 
t. The moments have been rescaled with U = 3 for visibility. 
Model B with a = 0.5, e d = 0, 7 = 0.5, T = 0 and 14 = 0.5. 


B. Magnetization 

In Fig. fl7j we plot the magnetization M z = ( Q^ — 
Ql)/ 2 as afunction of the magnetic field h for various 
strength of the interaction. At small field one expects 
M z oc Ii/Tk where the Kondo temperature Tk is the 
characteristic scale of the Kondo effect. We do observe 
indeed a sharp rise of the spin susceptibility when we 
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FIG. 15. First 6 moments of the sign s n associated with 
the currents I n versus time t. Model B with a = 0.5, td = 0, 
7 = 0.5, T = 0 and 14 = 0.5. Assymptotic s6 = 0.025. 



FIG. 16. First 6 moments of the weight c n associated with 
the currents I n versus time t. Model B with a = 0.5, td = 0, 
7 = 0.5, T = 0 and 14 = 0.5. 


switch on the electron-electron interaction. As a consis¬ 
tency check, the full lines in Fig.[l7]reproduce the results 
obtained in Ref. 133 by solving the Bethe-ansatz equa¬ 
tions. We find a remarkable agreement between the two 
(supposingly exact) techniques. Note that the agreement 
is not supposed to be perfect as the Bethe-ansatz equa¬ 
tions assume the universal regime r <C D (D = 4 is the 
band width) while the QMC calculations are performed 
for the microscopic model with a finite value T = 1/4. 
Note that, as in the previous section, the results have 
been calculated at zero temperature and bias voltage 
which is the most difficult situation. Indeed, at finite 
temperature the non-interacting Green’s function decays 



FIG. 17. Zero temperature magnetization. M z = (Q ^ — 
Q\)/2 as a function of magnetic field h for model A with 
td = —0.25, 7 = 1/4, T = 0 and t = 50. Lines are analytical 
results (from Bethe-ansatz Ref l37l) . in particular red line is 
at U = 0, black line at U = 0.7r (r = 47 2 = 1/4), blue 
line at U = 1.5F, dashed lines with points are results from 
QMC. The actual calculations were made with model B with 
a = 0.5 (h > 0.025) and a = 0.8 (h < 0.025), while enforcing 
t d = aU — 0.25. 

faster with time, which ensures a faster convergence of 
the results with time and consequently smaller values of 
c n . Following the strategy discussed in the previous sec¬ 
tion, the calculations are performed with a finite value of 

a. 


VIII. DISCUSSION 

In summary, we have presented a general purpose al¬ 
gorithm to calculate systematically the electron-electron 
interaction corrections to the physical observables of a 
generic nanoelectronic circuit. We validate the approach 
using the non-equilibrium Anderson model and find that 
we can calculate the first moments (up to n = 15 in this 
work using a recent laptop computer) in the stationary 
regime and without being plagued by the appearance of 
the dynamical sign problem. In a second step, one evalu¬ 
ates the interacting series whose radius of convergence is 
a priori unknown. Our results indicate that this radius 
of convergence is strongly affected by the mean-field part 
of the interaction (role of the a parameter). 

There are many aspects left to future work. Clearly, 
one wants to apply the technique to other, larger sys¬ 
tems such as a quantum point contact (”0.7” anomaly) 
or a quantum dot embedded in an interferometer. The 
technique could also be directly generalized to address 
several particle species (fermions or bosons), hence to 
study models coming from, e.g. circuit QED. Another 
route is to design techniques to obtain series with larger 
radius of convergence in order to reach more correlated 
regimes. This could be achieved using standard series 
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analysis tools (Borel resummations...) but also by com¬ 
puting the perturbative series for other quantities, such 
as the self-energy. Last, one needs to develop the mea¬ 
surement of more complex objects, such as the Green’s 
function, in order to be able to use the present technique 
for self-consistent schemes like DMFT or cluster DMFT. 
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first column, one finds det R„ = aidet M/. R„ and M„ 
having all their columns equal but the first one, they can 
now be put together into one determinant of a single ma¬ 
trix A = det M„ — det R„ = det [M ra + L„] (which is 
equal to M„ except for its upper corner shifted by — aq). 
This procedure can be continued with 0 . 2 , 0 .%... until no 
Oi are left which proves the above statement. Note that 
in the particular case considered in this work, the ma¬ 
trix M„ + L„ is block diagonal with respect to spin so 
that its determinant factorizes into two smaller determi¬ 
nants. This can be used for a faster calculation of the 
said determinants. 


Appendix B: Clustering property of the sum over 
Keldysh indices 


Appendix A: matrix M n for density-density 
interaction 

In this Appendix, we study the form of the matrix M„ 
(defined in Eq. ©) when the electron-electron interac¬ 
tion takes the form of a density-density interaction: 

Hint(f) = Vjj (f)[c|cj — Oi\[c^Cj — Oj] (Al) 

ijkl 

This form is more restrictive than the one studied in the 
main body of the manuscript. However, it allows for 
some optimizations as well as the inclusion of the one- 
body correction (proportional to cq) which have proved 
to be important and is sufficient for the applications to a 
single interacting site shown in this article. The result is 
very simple: the matrix M„ must be replaced by M„+L„ 
where the diagonal matrix L n consists of the cq: 


In this section, we show that summing over the 
Keldysh indices leads to a clustering property, i.e. that 
when the times of the integration in Eq. (14) are far 


from the time t where the observable is computed, the 
integrand decays. Let us consider Eq. ( fl5| ), in the case 
where e.g. u\,...,u p are far from t, e.g. close to t — At 
and u p +i,...,u n are close to t. We study the case 
where At becomes large. Let us denote u p +i,... ,u n ,t 
by v \,..., v n - v to simplify the notations. Then the ma¬ 
trix of Eq. (15) has a block structure of the form: 


M„ = 


A B 
C D 


(Bl) 


where 


■^i,j — 9 1 ^j ) 

Ci,j = g(vi,Uj) 


Dl,m = g(vi,v m ) (B 2 ) 


L — 7 . 



/ a h 

0 

0 

0 ... 0\ 


0 

a j 1 

0 

0 ... 0 


0 

0 

0^2 

0 ... 0 

l 

0 

0 

0 

Oj 2 ... 0 





. 0 


V 0 

0 

0 

0 ... 0 J 

from Ref. HI 

We present here 


(A2) 


recursive proof for completeness, 
averages of the form 


We need to evaluate 


A = (( e lci - ai)(c£c 2 - a 2 ) • • • (cjyfov - a N )) (A3) 


When all the cq vanish, the above average is equal to 
A = det M n . Let us now switch on u\. The previous 
result becomes A = det M„ — aidet M/ where the ma¬ 
trix M/ is identical to M„ with its first line and row 
removed. Equivalently, M/ can be replaced by a ma¬ 
trix R n of the same size as M„ with the first column 
filled with (ai, 0 ,..., 0 ) and all the other columns are 
those of M„: one immediatly checks that upon develop¬ 
ing the corresponding determinant with respect to the 


where i, j = 1,... ,p , and l,m = 1,... n — p. Therefore 

detM„(C n , {cq}) = det Adet^D — BA _1 C) (B3) 

Our assumption is that |iq — Vj\ ~ O(At) for At —» 00 . 
Moreover, at large time the non-interacting Green’s func¬ 
tion decays, as 1/t. Since B and C contain only Green’s 
functions with one u and one v, their arguments are of 
order At, and the matrix elements decay as 0(1/At). 
Hence the matrix M n is block-diagonal and 

det M„(C n , {oi}) = det A det D + 0(1/At 2 ) (B4) 

A is in fact a P p matrix, which leads to 

(-l) E? = lQl detA = 0 (B5) 

ai,...,a p 

Moreover D does not depend on the first p Keldysh in¬ 
dices, which gives finally 

Y, (-l) Er = ia ’detM„(C„,{aJ) = 0(1/At 2 ) (B6) 

a 1 . ,a n 
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We see from this argument that at large separation (for expansion to commute. 

Ui far from i), the terms at fixed Keldysh indices do not 

decay, while the sum over the Keldysh indices does. This 

is consistent with the observations made in Sect. hyi a 

, £ ,, rr- . , ^ , Appendix C: Properties of M\C n \ 

this is also necessary for the coefficients Q n to converge 1 1 

at long time, i.e. for the long time limit and the U- 


An important property of the M„ is found by considering the role of the ” mirror” configurations where the Keldysh 
indices {cii} are replaced by {1 — a,}. Let us consider the local density first: c|(t)cj(t). We suppose that the iq are 
ordered from the smallest to the largest (if not we relabel them). Then, using Wick theorem, we have 


det M n (C n , {a,}) = i 2 " +1 (H int (u 1 ) ajv H int (u 2 ) a2 ... H int (u„) ai cj(t)c i (t)H int ('U„) 1_OJV 
from which we immediately find that 


■ ■ ■ H int (u 2 ) 1_a2 H int (u 1 ) 1_ai ) 

(Cl) 


det M„(C„, {a,}) = - det M„(C„, {1 - a*})* 


(C2) 


By adding together a configuration and its mirror, we deduce from the above that 

- i n+1 Y, °* det M ^Cn, {*}) 

{a;} 


(C3) 


is real and therefore M[C n ] from Eq. ( f29| ) is real as well, 
replaces c\(t)ci(t) by i[c}(t)cj(t) - ct(t)c»(t)]. 


Appendix D: Gray code for fast updates of the 
determinants 

The sum over the Keldysh indices can be accelerated 
using a Gray code, see e.g. Ref. [38] section 20.2. Using 
a Gray code, it is possible to enumerate all the integers 
from 0 to 2 N — 1 in such a way that their binary repre¬ 
sentation change only by one bit at each step, and that 
the last one contains only one bit (hence is also one bit 
flip away from 0). To do this, we start with 0, and for 
n = 0, ..2 N — 1 we flip the bit at position c where c is 
given by: 


j ffs(~ n) iif n < 2 N — 1 
“={« iif„ = 2»-l (D1) 

where the ffs(z) returns the position of the first (least 
significant) bit set in the integer i. and ~ n is the binary 
complement of the integer n (using the syntax of the C 
language). We obtain all integers between 0 and 2^ — 1 
once and only once. The last case ensures that the integer 
returns to 0. 

The Keldysh indices at order N are a list of size N 
of 0 (upper contour) and 1 (lower contour). They are 
therefore in one-to-one correspondence with the integers 
from 0 to 2 n — 1, via their binary representation. Using 
the Gray code, we can enumerate the Keldysh indices 
by changing only one bit at a time, which in our algo¬ 
rithm means changing one line and one column in the 
determinant at a time. This operation is of complexity 
iV 2 , i.e. much quicker than a full recomputation of the 


The same argument holds for the current provided one 


determinant N 3 , and it is very commonly used in all de- 
terminantal Quantum Monte-Carlo. It is implemented 
using BLAS 2 operations. 

The following piece of C++ code illustrates the use of 
the Gray code in our implementation: 

auto two_to_N = uint64_t(l) << N; 

for (uint64_t n = 0; n < two_to_N; ++n) { 
int nlc = (n < two_to_N - 1 ? ffs(~n) : N); 

// Change the line and column numbered nlc 
// in the matrix 

> 

where N is the the order at which we compute the de¬ 
terminant, Note that the matrix has returned after the 
loop to the value it had before the loop, up to numerical 
errors, which are typically controlled at this place. 


Appendix E: Extrapolating series beyond their 
convergence radius 


In this appendix, we use the Lindelof 32 ^ extrapolation 
method, for the computation of the charge Q{U) with 
a = 0. On Fig. 0 we have seen that for a = 0, the 
series can not be summed for U > 0.6 directly. Using the 
Lindelof formula: 


N 


Q(U,N , e) = YQnU n e~ mln ^ 


(El) 


n —0 


we plot Q(U, N, e) for U = 1.5 as a function of e in Fig. 18 
We see that the extrapolation method produces a result 
for Q{U = 1) in very good agreement with the imaginary- 
time QMC (dashed line in Fig. Ill result. 






19 



FIG. 18. Extrapolation of the series. Q{U = l,e) as a 
function of e for model A with = 0, T = 0 and 7 = 0.5 
for various maximum order N = 4,5, 6 , 7. The circles show 
a simple extrapolation from a linear regression in the region 
where the various curves overlap. We do recover Q(U = 1) « 
0.32 well beyond the apparent convergence radius of the series. 


Appendix F: Role of bound states in model A 


which exists when || > 2(1 — 7 1 2 3 ). The top panel of 
Fig. [19] shows the corresponding bound state energy cal¬ 
culated numerically (which matches perfectly the above 
expression) as a function of e^. When 7 < 1, there is a 
finite window of values of €d where there are no bound 
states in the system. The bottom panel of Fig. [l9| shows 
the associated non-interacting lesser Green’s function of 
the system as a function of time. In the absence of bound 
state (dashed line, 7 = 0.5), we see that gQ 0 {t) decays to¬ 
wards zero at large time. Note that this decay would 
be much faster a larger temperature. In the presence of 
bound states however (full line, 7 = 1 ), we find that at 
large time Qq 0 saturates to its bound state contribution 
gQ 0 (t ) oc e~ lElt . As a results, the convergence of the 
terms in the perturbative expansion with the time t will 
be much slower and are not even guaranteed to converge 
as the bound states do not relax. This is illustrated in 
Fig. [20] which shows Q 2 (defined in Eq. (28 1 ) versus 1/t 
for 7 = 1 (upper panel, a bound state is present for any 
ed 7 ^ 0) and 7 = 0.5 (lower panel, no bound state). Sim¬ 
ilarly, Fig. [ 2 l| shows Q 2 versus €d with (lower panel) and 
without (upper panel) bound state. We find that without 
bound state, Q 2 quickly converges towards its stationary 
value (roughly for t > 10T). However, in presence of 
bound states the convergence is much slower. 


As an illustration of the role of bound states, we briefly 
study their role in model A. For this model, they can be 
obtained analytically by simple wave matching: One gets 


E 1 


2q 2 - 1 


(7 2 


l)kd| + l 2 \J^d + 4 * 6 * 8 (2t 2 — 1) 
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